Numerical simulations of the decay of primordial magnetic turbulence 
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We perform direct numerical simulations of forced and freely decaying 3D magnetohydrodynamic 
turbulence in order to model magnetic field evolution during cosmological phase transitions in the 
early Universe. Our approach assumes the existence of a magnetic field generated either by a process 
during inflation or shortly thereafter, or by bubble collisions during a phase transition. We show 
that the final configuration of the magnetic field depends on the initial conditions, while the velocity 
field is nearly independent of initial conditions. 

PACS numbers: 98.70.Vc, 98.80.-k 



I. INTRODUCTION 

Astronomical observations show that galaxies have 
magnetic fields with a component that is coherent over a 
large fraction of the galaxy with field strengths of order 
10~ 6 Gauss (G); see Refs. [l|-[4| and references therein. 1 
Modeling the origin of these fields is a challenging pro- 
blem. Here we consider models in which the seed field is 
generated in the early Universe, see, e.g. 0, H, @. There 
are a number of such models and different models result 
in different magnetic field structures. 

It is well known that quantum-mechanical fluctua- 
tions during inflation can generate energy-density inho- 
mogeneities that seed observed large-scale structure, see, 
e.g., Ref. 0. Quantum- mechanical zero-point fluctua- 
tions can also generate a seed magnetic field, provided the 
relevant abelian gauge field is such that the lagrangian 
density is not conformally invariant during inflation @- 
[Toj . A seed magnetic field generated during inflation 
has a correlation length today that can be very much 
larger than the current Hubble radius. A convenient phe- 
nomenological way of breaking conformal invariance dur- 
ing inflation is to couple the inflaton scalar field <f> to the 
vector abelian gauge field by e a '^F ia ,F ,iV , where F^ is 
the vector field strength tensor and a is a parameter [9j . 
Depending on the value of a, this can result in a suffi- 
ciently large seed magnetic field to explain the observed 
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1 See Ref. and references therein for a technique that might 
soon prove useful for measuring a larger-scale, cosmological, mag- 
netic field. 



galactic magnetic fields. This is a classically consistent 
model, with vector field back-reaction during inflation 
being negligibly small 0, |Tl| . While the abelian gauge 
field becomes strongly coupled at early times [tj dH [l2| , 
this is not a problem for the effective, phenomenological, 
classical model, @, [ill an d also see d, HH . Of course, 
just as one does not reject the effective, phenomenologi- 
cal, classical "standard" ACDM cosmological model, one 
does not reject inflation-based seed magnetic field gener- 
ation models for being quantum-mechanically inconsis- 
tent. Rather, much as the case for ACDM, it is of great 
interest to understand how such a successful inflation- 
based seed magnetic field generation model might arise 
from a more fundamental underlying model or theory. 

There are a number of other seed field generation mod- 
els; see Ref. [l4j for a discussion of seed field generation 
in Hofava-Lifshitz gravity. Yet another interesting pos- 
sibility is primordial magnetic field generation during a 
cosmological phase transition [15l ] . If the phase transition 
takes place at late times the magnetic field correlation 
length is smaller than the Hubble radius. 

If the magnetic field was generated during or shortly 
after inflation (during a pre- or re-heating phase transi- 
tion), and survived to the epoch of recombination (the 
last-scattering surface), it should leave observable signa- 
tures on the cosmic microwave background (CMB) ra- 
diation fluctuations, see Ref. [l6j] for a review and Refs. 
[13, EH f° r more recent discussions. 2 

The shape and magnitude of a primordial magnetic 
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field can be constrained observationally. Magnetic field 
energy density scales like radiation (i.e. relativistic) en- 
ergy density. 3 During nucleosynthesis, any new form of 
radiation-like energy density is constrained observation- 
ally to be less than about 10% of the usual radiation en- 
ergy density. Hence, agreement between the big bang nu- 
cleosythesis (BBN) model light nuclei abundance predic- 
tions and the observations yields a constraint on the pri- 
mordial magnetic field energy density: ris^o — 2.4xl0~ 6 
[2l| . Here Qb is the current value of the magnetic field 
energy density parameter and ho is the current value of 
the Hubble parameter in units of 100 km s _1 Mpc -1 . 
Naturally, this limit holds only if the primordial magnetic 
field was generated prior to or during nucleosynthesis. 

Other constraints arise through the effects of a pri- 
mordial magnetic field on the propagation of CMB pho- 
tons (Faraday rotation of the CMB polarization [l8| , and 
magnetic field induced scalar, vector, and tensor modes 
of CMB anisotropies jla. Il7l|). Available data limit the 
current value of a cosmological magnetic field to be less 
than about 10~ 8 G for a scale-invariant or homogeneous 
primordial magnetic field. 

Another interesting signature of a primordial mag- 
netic field is relic gravitational waves generated by the 
anisotropic magnetic stress [22]. The amplitude of the 
gravitational wave signal is determined by the efficiency 
of gravitational wave generation and the strength of the 
primordial magnetic field. The efficiency of gravitational 
wave production is small [23| due to the small value of the 
Newton constant G in the coupling between the primor- 
dial magnetic field and the gravitational waves. However, 
a magnetic field that satisfies the BBN limit and that is 
generated during a strong-enough first-order phase tran- 
sition can lead to a detectable signal for LISA. Therefore, 
the direct measurement of the gravitational wave back- 
ground can lead to an independent limit on a magnetic 
field generated in the early Universe [24j . 

On large scales, because of the high conductivity of 
the plasma in the early Universe, the magnetic field 
is treated as a frozen-in field with its evolution deter- 
mined by the simple dilution of magnetic field lines, 
B(x, t) oc B(x)/a 2 (i), where t is the physical cosmic time 
and a(t) is the scale factor. In general, however, the evo- 
lution of a primordial magnetic field is a complex process 
influenced by MHD as well as by cosmological dynam- 
ics [25[. In particular, the presence of a magnetic field 
can dramatically affect the behavior of primordial tur- 
bulence (for example, turbulence associated with phase 
transition bubble motions) |26l - l28l |. Also, the presence of 
a magnetic field might itself lead to the development of 
turbulent motions, see Refs. (29l - [3l| and Sec. Ill below. 

In this paper we study the evolution of a primordial 
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magnetic field that is coupled, through the MHD equa- 
tions, with the fluid, during a cosmological phase tran- 
sition. We consider two different possibilities: (i) when 
the magnetic field energy has been injected in the plasma 
with no initial vorticity perturbations (no turbulent mo- 
tions) at a given typical scale; and, (ii) when vorticity 
perturbations are present during the phase transition and 
couple to the magnetic field with a given spectrum. Both 
these assumptions can be justified. The magnetic field 
might be generated prior to the phase transition, with 
the generation mechanism not requiring turbulent mo- 
tion. On the other hand, the presence of initial turbulent 
motions can also be justified since bubble nucleation, ex- 
pansion, and collisions can lead to primordial kinetic tur- 
bulence [32l |. Of course, in reality both magnetic field and 
kinetic turbulence generation can take place together, 
while in our approach we assume the a priori presence of 
magnetic energy and/or kinetic energy during the phase 
transition. 

A causal, phenomenological description of MHD on 
large scales during the cosmological phase transition mo- 
tivates a white noise spectrum for the magnetic spectral 
energy density, £7m(&) oc k 2 , where k is the wavenumbcr 
[33d . A steeper Batchelor spectrum with EM{k) oc k 4 has 
been claimed to follow from cosmological causality and 
the divergenceless condition of the magnetic field [34j . 
To understand the evolution of a cosmological magnetic 
field, and to be able to make observable predictions in 
this model, it is important to resolve this impasse. One 
would like to know whether the final large-scale magnetic 
energy spectrum evolves to something steeper than white 
noise, or to a spectrum closer to the Batchelor one. 

In fact, the MHD process itself might establish a dif- 
ferent spectrum: The possibility of generating a random 
magnetic field from isotropic turbulence was first pro- 
posed by Batchelor (35|. He invoked the imperfect anal- 
ogy between vorticity and magnetic fields that should 
generally imply a large scale distribution of the magnetic 
energy similar to the kinetic energy of turbulent motions 
Ex(k) oc k 4 (36|. Later, Kazantsev [37j was able to rig- 
orously establish the possibility of small-scale dynamo 
action. He assumed that the velocity field varied only 
on large scales and found that weak magnetic fields are 
amplified mainly on the resistive scale. Initially, this was 
thought to be applicable only to turbulence at large mag- 
netic Prandtl number [38| , where the viscous cutoff scale 
is much larger than the resistive cutoff scale. Kazantsev 
found that the magnetic energy spectrum increases with 
wavenumber like fc 3 / 2 , which is slightly shallower than the 
white noise spectrum k 2 . The emergence of a Kazantsev 
spectrum turned out to be much more ubiquitous and not 
only applicable at large magnetic Prandtl number. Sim- 
ulations at magnetic Prandtl numbers of order unity also 
clearly showed the Kazantsev k 3 ^ 2 spectrum 27]. Such 
a spectrum could be of interest for primordial magnetic 
field evolution, because it implies somewhat larger power 
at large scales than white noise. 

The main goal of our paper is to determine, through 
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MHD modeling, the evolution of the magnetic energy 
spectrum. To keep our study as general as possible we 
do not make any assumption about the physical process 
leading to primordial magnetic field generation. Also, 
we do not address the phase transition physics itself, 
keeping the total magnetic field energy density ps as a 
free parameter whose maximal value is fixed by the BBN 
bound. Obviously, if the magnetic field is generated dur- 
ing a phase transition, pb will depend sensitively on the 
amount of latent heat that is transformed to magnetic 
energy, i.e. on the efficiency of the magnetic field gen- 
eration process. One important issue addressed here is 
to determine the spectral shape of the magnetic field at 
large scales, assuming that magnetic energy and vortic- 
ity perturbations are closely coupled during the phase 
transition. Another important question is related to the 
duration of MHD turbulence and how long it takes to 
reach equipartition between kinetic and magnetic energy 
densities in the primordial plasma. 

The structure of this paper is as follows. In Sec. II we 
describe the model, defining the main model parameters 
(Sec. HA), characterizing the magnetic field spectrum 
(Sec. IIB), and formulating initial conditions (Sec. IIC). 
In Sec. Ill we present results from direct 3D MHD simu- 
lations performed using the Pencil Code [39j | . We con- 
clude in Sec. IV. We employ natural units with h = 1 = c 
and gaussian units for electromagnetic quantities. To 
properly account for the expansion of the Universe we 
use comoving quantities with conformal time t that is 

1 /2 

related to physical time t p h ys as t oc t' h during the ra- 
diation dominated epoch. 

II. MODEL DESCRIPTION 

We assume that magnetic energy is generated at the 
electroweak or QCD phase transition, or during inflation. 
In the first case magnetic energy is explicitly injected into 
the fluid on small length scales, smaller than the Hubble 
length at the moment of the phase transition. If the mag- 
netic field originated at an earlier epoch of inflation the 
length scale at which the magnetic field interacts with 
the fluid is set by the characteristic length scale of the 
system, and is again smaller than the Hubble radius. In 
the absence of magnetic or kinetic helicity we do not ex- 
pect an inverse cascade (i.e., energy flow from smaller to 
larger length scales). 

A. Phase transition characteristics 

To model magnetic field evolution one needs to know 
the physical conditions during the phase transition. 
First, generation of turbulence requires a first-order 
phase transition so that at the critical temperature of 
the phase transition bubbles of the new vacuum nucleate 
within the false vacuum. Bubble collisions then gener- 
ate turbulent motions. The standard electroweak model 



does not have a first-order phase transition , and can- 
not account for baryogenesis. However, modifications of 
the standard model, such as the minimal supersymmetric 
standard model (MSSM) [4l| . result in first-order elec- 
trowe ak p hase transitions and can account for baryoge- 
nesis [421 ]. Also, recent lattice QCD computations [43[ 
have not yet excluded that the QCD phase transition is 
a first-order one, with bubble nucleation and collisions. 

The main parameter characterizing turbulent motions 
is the r.m.s. velocity vq which determines the kinetic en- 
ergy density of the turbulence. Obviously, vq depends 
sensitively on the phase transition physics, and in partic- 
ular on the phase transition bubble wall expansion veloc- 
ity Ufc, [44]. To model the development of turbulent mo- 
tions during the phase transition we adapt earlier analyt- 
ical or semi-analytical results 0, IH| . The first question 
that must be answered is whether the phase transition 
is first order, and, if so, what fraction of total available 
vacuum energy is transformed into kinetic energy of the 
bubbles. The r.m.s. velocity of the turbulent motions can 
be approximated as (see Ref. [45j]) 

Here a — /0 V ac/pthcrmai is the ratio of the vacuum energy 
density associated with the phase transition to the ther- 
mal energy density of the Universe at the time (a char- 
acterizes the strength of the phase transition), and the 
efficiency factor n is the fraction of the available vacuum 
energy that goes into the kinetic energy of the expanding 
bubble walls (as opposed to thermal energy). 

On the other hand, the primordial magnetic field is 
characterized by the r.m.s. Alfven velocity given by 

va = /A B = X [i^ ^ 7.65 x 10 2 v/fi^. (2) 

Here w rad = 4p rad /3 (with p rad ~ pthcrmai) is the radi- 
ation enthalpy of the relativistic fluid and we have used 
^rad/i 2 ] = 2.56 x 10 -5 for a present-day CMB temperature 
T = 2.74 K. At temperature T+ of the phase transition 
we have p r ad(T*) = 7r 2 (/*(T*) /30, where g± is the num- 
ber of relativistic degrees of freedom at temperature 
The r.m.s. Alfven velocity does not depend on T*, but it 
is weakly dependent on g+, 

Assuming equipartition between kinetic and magnetic en- 
ergy densities, one has va — v$. While the equipartition 
condition can be justified by MHD dimensional analy- 
sis [29j, we explicitly show that it holds by performing 
a 3D direct numerical simulation of a primordial mag- 
netic field coupled to fluid motions; see Sec. III. Also, 
partial equipartition, vq ~ 0.8va, is reached in the case 
where the magnetic field was generated first at the in- 
jection scale, which then led to rapid growth of vorticity 
perturbations in the initially no-turbulent plasma. 
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The bubble wall expansion sets the maximal size of the 
bubble, which we associate with the size of the largest 
turbulent eddy, l = Vb/3 -1 , where (3 is a parameter that 
characterizes the duration of the phase transition. In 
particular, /3 can be obtained from the bubble nucleation 
rate (32J. The Hubble time H^ 1 at the phase transition 
is another characteristic time, and it sets the "causal- 
ity" horizon. At this point it is useful to define the pa- 
rameter 7[= loH+ = Vb(f3/H ± ")~ l \ which determines how 
many maximal-sized bubbles are within the Hubble ra- 
dius, N w 7 3 . 




FIG. 1: The maximal allowed value _Em(&) of the magnetic 
field generated during the EW phase transition with T* = 100 
Gev, = 100, and 7 = 0.01. 



B. Magnetic field spectrum 



If a primordial magnetic field is randomly oriented and 
its mean value vanishes, i.e. (B(x)) = 0, it is conveniently 
described statistically in terms of the n-point field corre- 
lation functions. If the field is isotropic with a Gaussian 
distribution, the magnetic field characteristics are com- 
pletely determined by the two-point correlation function 
(B(x + r, t + t)B(x, t)). To construct this main charac- 
teristic function we need to know the spatial distribution 
(i.e., the correlation length) and the temporal evolution 
of the magnetic field. 

When considering a causally generated primordial 
magnetic field, its maximal co-moving correlation length 
Cmax is set by the co- moving Hubble radius Xh(= 
H^clq/ a*) at the moment of generation (here 00 and 
a* are the scale factors today and at magnetic field gen- 
eration, respectively). Causality implies that 



lation length is to use the magnetic energy spectrum 4 
EM(k,t), defining the correlation length by £m(£) = 
[Jo 00 dk k- l E M {k, t)] I r °° dk E M (k, t). 

In most models of magnetic field generation during 
a phase transition (To| . the magnetic field correlation 
length is determined by the phase transition's bubble 
sizes. In this case a characteristic magnetic field corre- 
lation length is assumed to be determined by the largest 
bubble, Ao — 7A#. This simple dimensional descrip- 
tion implies that the magnetic energy density is redis- 
tributed through MHD processes and establishes a mag- 
netic spectrum with spectral energy density measured 
today that is given by [23j 



E M {k) < 
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Here, k = k/ko and kp = kjj/ko with fco = 2t:/Xq, 
and kjj = 2tt/\]j is the damping wavenumber deter- 
mined through the viscosity-driven dissipation of the 
magnetic field. In the case of stationary non-helical 
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(4) 

where the temperature corresponds to the energy scale 
at field generation. For the electroweak phase tran- 
sition, T* is related to the Higgs mass Mh through 
T* ~ (1.2 ± 0.2)Mh, so we parameterize the depen- 
dence by normalizing to a temperature of 100 GeV. Also, 
for the electroweak phase transition g± ~ 100, but the de- 
pendence on g+ is much weaker than the T± dependence. 

Another way to determine the magnetic field corre- 



4 In what follows we use [46| 

(B*(k,t)B,(k',t + r))=J(k-k')^(k,t)/[#),r], (5) 



where 



87rfc 2 F"(k,r) 



: 2P !J (k)E M (M) + is i:i ikiH M (k,i) . 



(6) 



Here Pij(k) = <5,;j —kikj/k 2 is the projection operator, Sij is the 
Kroncckcr delta, k = |k|, Eiji is the totally antisymmetric ten- 
sor, and r)(k) is an autocorrelation function that determines the 
characteristic function f[rj{k), t] describing the temporal decorre- 
lation of turbulent fluctuations. The scalar function HM(k,t) is 
the magnetic helicity spectrum. All configurations of the helical 
magnetic field must satisfy the "realizability condition" 29, 30]. 
\f™dkH M (k,t)\ < 2f M (t)£ M (t). 
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Kolmogorov turbulence the Reynolds number determines 
the damping scale ko as Re = (k D /k ) 4 / 3 . Note that 
the Reynolds number is high enough in the early Uni- 
verse to ensure the presence of a wide inertial (turbu- 
lent) range, ko < k < fc^. 5 The large-scale behavior 
of the magnetic field is determined by the parameter a. 
The scale- invariant spectrum corresponds to a = — 1 Q, 
the Kazantsev spectrum has a = 3/2 [37j . the white noise 
spectrum corresponds to a = 2 [33| , and the steep Batch- 
elor spectrum has a = 4 [35]]. Equation ([7]) does not ac- 



count for any damping of magnetic energy and can be 
viewed as the BBN bound imposed at the moment of the 
establishment of the magnetic spectrum; see Fig. 1. 



C. MHD formalism and initial conditions 

As discussed in Ref. [13], the usual relativistic MHD 
equations are identical to the MHD equations in an ex- 
panding Universe with zero spatial curvature when phys- 
ical quantities are replaced by their co-moving counter- 
parts and conformal time rj is used in place of physical 
time. Based on this fact, we perform direct numerical 
simulations of MHD turbulence in the expanding Uni- 
verse using the usual MHD equations for an ultrarela- 
tivistic equation of state, but nonrelativistic bulk mo- 
tions. 

The characteristic length-scale of the problem is set by 
the phase transition bubble size 'fXn ■ The typical time- 
scale is the eddy turn-over time. Assuming that turbu- 
lent eddies correspond to phase transition bubbles, the 
physical turn-over time is tq ~ (vb/vo)(/3/-ff*) -1 , where 
«0 = (V 2 ) 1/2 - 

Another characteristic of the initial stage is the amount 
of magnetic and kinetic energies present during the phase 
transition, i.e. the initial Alfven and turbulent r.m.s. ve- 
locities. As discussed above, both quantities are sensi- 
tively dependent on the available vacuum energy that is 
converted to magnetic energy (if the magnetic field was 
generated during the phase transition) and/or the kinetic 
energy of the turbulent motions. In our analysis below 
we assume that about 10% of the vacuum energy is in 
the form of initial magnetic energy, which corresponds to 

VA,in ^ 0.3. 

We perform the numerical simulations of magnetic field 
evolution in two stages. During the first stage, we model 
forced MHD turbulence with injection of energy at fixed 
wavenumber. This type of driving force is supposed to 
mimic the action of bubble-induced external forces in 
MHD turbulence during the phase transition. The simu- 
lation of forced MHD turbulence is carried out before 
equipartition is reached, i.e., before vq ~ Va- After 
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number in the early Universe at the scale of energy injection of 
the turbulence and magnetic field. 



equipartition is reached we switch off the driving force 
and allow free decay of the turbulent state. 

We use two different types of external forcing during 
the first stage of the simulations. These are injecting 
magnetic or kinetic energy at a given scale associated 
with the turbulent eddy size. These cases differ in the 
initial magnetic field configuration. Injection of the mag- 
netic energy in the flow is achieved by using a delta func- 
tion spectral energy density function for the magnetic 
field. Injection of the kinetic energy is achieved by us- 
ing a delta function spectral energy density function for 
the velocity field. In the latter case the spectral energy 
distribution of the initial magnetic field is of Batchelor's 
type (Eb oc k 4 ) and its amplitude is close to that of the 
kinetic energy density (the equipartition condition). 

We also define the spectrum of the velocity field 
Eic{k,t) and the total energy density £ic(f) — 
f °° E K (k,t)dk. One of our main goals is to determine 
whether the presence of magnetic fields in the early Uni- 
verse (e.g. generated prior to a phase transition) can lead 
to strong turbulence. 



III. NUMERICAL SIMULATIONS 

Numerical simulations of the magnetic field evolution 
were performed using the Pencil Code [27 l [48l - 450l | ; see 
[39j | for the website. 

We perform all simulations using co-moving quantities 
and conformal time. For simplicity we work with dimen- 
sionless quantities, and use k\ — fco/30 as our wavenum- 
ber unit. The chosen box size covers a wavenumber 
range from fco/30 up to 4.3 x ko, and the maximal length 
scale considered corresponds to Ai = 30-f\n, which is 
still within the electroweak phase transition Hubble scale 
(7ew < 0.01). In the case of a QCD phase transition for 
an extremely large QCD bubble velocity, Vf, — > 1, the 
scale can exceed the QCD phase transition Hubble scale 
by a factor of 2 (7qcd < 0.15). Of course, in the sim- 
ulations we have to use a relatively large value for the 
dissipation wavenumber, which is fcrj ~ 2ko at the end 
of the simulation. This high value is a consequence of 
choosing a constant viscosity that must be large enough 
so that it can also cope with the initially much larger 
value of the energy dissipation rate. Such high values are 
obviously not realistic for the early Universe where the 
Reynolds number is extremely high. However, we moti- 
vate this choice by the fact that we are mainly interested 
in the evolution of the magnetic field outside the inertial 
range for k < ko- The time unit in our simulations is 
set by the computational box size and sound speed, i.e., 
h = y/3/ki. 

As is well known, free decay of MHD turbulence im- 
plies an increase of the magnetic eddy size with decreas- 
ing magnetic energy density. Figure [2] illustrates this 
fact in our simulations. We display the y component of 
the magnetic field on the periphery of the domain dur- 
ing stages of its evolution after equipartition has been 
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reached and driving was switched off, in the case when 
the magnetic energy was injected at some typical scale. 

The characteristic time scales for the qualitative 
changes in the magnetic field distribution is approxi- 
mately equal to 20ii ~ 1657A// - a value that exceeds 
slightly the Hubble time-scale, Xh/c, for the electroweak 
phase transition. 



A. Development toward equipartition 

We first examine how long it takes for forced MHD 
turbulence to establish equipartition between the mag- 
netic and kinetic energy densities, i.e. where va — Vq. Of 
course, the evolution of the magnetic energy spectrum is 
strongly scale dependent. We first address the case when 
the magnetic field is "created" at a typical scale k^ 1 , 
corresponding to the largest bubble. The Alfven veloc- 
ity of this field is close to the maximal value allowed by 
the BBN bound, i.e. VA.in — 0.3. In this case 10% of 
the available vacuum energy has been somehow trans- 
formed into magnetic energy before or during the phase 
transition. We perform our simulations with zero ini- 
tial velocity perturbations. Such initial conditions may 
apply to the generation of turbulent motions during a 
phase transition via MHD processes, while the magnetic 
field was generated prior to that, for example through 
quantum fluctuations 0, [lOj . In this case the evolution 
of magnetic and kinetic energy spectra during the first 
stages of MHD coupling are shown in Fig. [3J Taking into 
account that the largest bubble size must be a typical 
length scale for turbulence, our initial conditions imply 
that at k = ho equipartition is reached almost instanta- 
neously. Larger scales need substantially longer times to 
establish equipartition. The initial evolution of magnetic 
energy spectra is shown in Fig. |4j 

In our phenomcnological description it has been as- 
sumed [45[ [5l[ that bubble collisions and nucleation lead 
to the development of turbulent motions, and the kinetic 
energy spectrum has been approximated as Kolmogorov- 
like in the inertial range (Ek oc k~ 5 ^ 3 ) and by a white 
noise spectrum (Ek oc k 2 ) for k < ko- As noted above, 
our simulations cannot adequately describe the inertial 
range and thus we cannot expect to see a Kolmogorov- 
like spectrum. On the other hand, the large scale con- 
figuration of the velocity field is well approximated by 
white noise; see Figs. 01 03 and [6] This agrees with what 
is predicted in the phenomenological approach. 

We note that the kinetic energy spectrum grows faster 
right after the driving is switched on. In particular, after 
ti, the r.m.s. velocity increases from zero to 0.17 while 
approaching 0.21 after 20ii. At vo ~ Q.8va we assume 
that equipartition is reached and, by switching off the 
driving, we allow the turbulence to enter the second, free- 
decay stage. 



B. Free decay of MHD turbulence 

The free decay of turbulence is shown in Figs. [5] and 
[5] for two different initial configurations of the magnetic 
field. In both cases we perform MHD simulations af- 
ter the driving force was switched off and free decay oc- 
curred. 

The first case corresponds to the process described 
above, so it was preceded by magnetic field injection into 
the turbulent plasma at the scale of the largest bubbles 
or turbulent eddies. The field was then allowed to de- 
cay, leading to the development of near equipartition be- 
tween the magnetic and turbulent energies. The second 
case corresponds to an initial configuration of a random 
magnetic field with a k A spectrum. In contrast to the 
first case, the process here can be roughly described as 
turbulent kinetic energy injection into the magnetized 
plasma. In this case we approximate the magnetic field 
spectrum at large scales by the steep Batchelor spec- 
trum (Em oc k 4 ). Of course, the real situation during 
the phase transition is somewhere in between. If the 
magnetic field was generated during the phase transition 
through bubble collisions [IH the same process of bub- 
ble collisions leads to the generation of turbulent mo- 
tions (vorticity perturbations). Thus, strictly speaking, 
we cannot split the turbulent motions and magnetic field 
generation and evolution. If the magnetic field was cre- 
ated before the phase transition, it affects the bubble 
collisions (locally inserting a preferred direction). As a 
result, the process of generation of turbulent motions is 
affected and this backreacts onto the magnetic field con- 
figuration itself. 

The first case — injection of magnetic energy into the 
turbulent plasma — provides a suitable setting for mag- 
netic field generation through the mechanism described 
in Ref. [52J , where the correlation length of the magnetic 
field is naturally set by the size of the phase transition 
bubbles. The second case with an established magnetic 
field spectrum is more appropriate for the causal mecha- 
nisms discussed in Ref. [34J. 

Figures [5] and show a distinctive difference between 
the turbulent states developed after the injection of the 
magnetic energy through a single mode magnetic field, 
and the injection of kinetic energy in the existing mag- 
netic field with smooth magnetic spectral energy density 
distribution. A major difference is revealed in the spec- 
tral distribution of magnetic and velocity fields at large 
scales. At large scales, the spectral slopes are approxi- 
mately 2 for kinetic energy and around 3 for magnetic 
energy. Therewith, the spectral distribution of kinetic 
energy shows less sensitivity to initial conditions and de- 
velops a shape close to white noise with Ek oc k 2 . The fi- 
nal configuration of the large-scale magnetic field slightly 
differs as a consequence of the initial conditions. In the 
case of magnetic energy injection (Fig. [5]) the spectral 
slope is shallower than 3, and tends to establish a white 
noise spectrum, while the case of kinetic energy injec- 
tion (Fig. [3]) most probably results in the steep Batchelor 



FIG. 2: Evolution of the turbulent magnetic field after turning off the forcing at time t = 14 1\. The B y component is shown 
on the periphery of the computational domain. 
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FIG. 3: Evolution of the magnetic energy spectral density 
EM(k,t) (upper panel) and kinetic energy spectral density 
Ei({k, t) (lower panel) for different values of k/ko = 0.03, 
0.06, 0.17, 0.33, and 1. The spectra are normalized such that 
jE M (k,t)dk = <u 2 )/2. 



spectrum. 

In all cases, throughout the free decay stage, the peak 
of the magnetic field spectral energy density drifts to 
smaller wavenumbers as fc pea k(t) ^ t~ l l 2 . At the same 
time magnetic power decreases as EM{k,t) oc t . Ac- 
counting for these scalings and noting that we use confor- 
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FIG. 4: The initial temporal evolution of magnetic energy 
spectra EM(k) [solid lines at t/ti = 0.1, 0.2, 0.5, 1.5, 5, and 
25, with (smoothed) Ej^(k) at k = lOfci increasing as t in- 
creases] are shown before the field reaches equipartition with 
the kinetic energy density. For comparison, kinetic energy 
spectra are shown for the same times (dashed lines). Thick 
lines (solid and dashed) indicate the last time, 25 1%. Straight 
dash-dotted lines have slopes 2 and 4. The box turnover time 
is 15 ti and is reached only for the last time shown. 



mal time, the temporal scaling in the case of an expand- 
ing Universe is somewhat slower than that in the case of 
laboratory grid turbulence where Ek oc t~£ with expo- 
nents between n = 1.13 [53| and 1.25 [Hij]. 

Simulations show that equipartition between kinetic 
and magnetic energies is sustained throughout the free 
decay stage of the turbulence. Therewith, at any given 
point in time, the kinetic energy spectrum peaks at the 
same wave-number as the magnetic energy spectrum. 
However, it seems that properties of the turbulent states 
developed through free decay depend on the method of 
their generation. The power spectrum of the turbu- 
lence developed after injection of a single-mode mag- 
netic field energy peaks at /c max ~ 3/ci, while injection 
of kinetic energy into an existing magnetic configuration 
leads to a turbulent state with power spectrum peaking 
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FIG. 5: Magnetic (solid) and kinetic (dashed) energy spectra 
in 12 regular time intervals of 4ii after having turned off 
the forcing, with (smoothed) spectra at k = 50fci decreasing 
as t increases, v — r\ = 10~ 4 in units of (kfti) -1 . The 
straight lines have slopes 3, 2, —2, and —1/2, with the first 
two near k = ki and the last two near k = lOfci. Thickest 
lines (solid and dashed) indicate the last time, which is 44 1\ 
since turning off the forcing. The intermediate thickness solid 
line, the highest or almost highest line for k/ki > 10, is the 
initial magnetic spectrum for this computation. 
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FIG. 6: Same as Fig. [5] but for a case where the initial mag- 
netic field had a k 4 spectrum close to equipartition with the 
velocity field, and then the forcing was turned off. Results 
are shown for nine times at intervals of 6ti. v = r\ = 10~ 4 
in units of (fcjti) -1 . The straight lines have slopes 2 and 3. 
Thickest lines (solid and dashed) indicate the last time, which 
is 48 1\ since turning off the forcing. The intermediate thick- 
ness solid line, the highest solid line for 5 < k/ki < 10, is the 
initial magnetic spectrum for this computation. 



at fcmax ~ 9fci. Hence, it seems that, in general, the char- 
acteristic length scale of the turbulence is sensitive to the 
initial driver: kinetic drivers result in smaller-scale tur- 
bulence states as compared to the case of magnetic field 
drivers. 



IV. CONCLUSIONS 

In this paper we have presented results from direct 
numerical simulations of primordial magnetic field evo- 
lution during cosmological phase transitions. These sim- 
ulations account for the expansion of the Universe. Sim- 
ulations were performed to model two different stages of 
primordial magnetic field evolution, the first stage when 
phase transition processes drive turbulence, and the sec- 
ond stage when free decay occurs. 

We show that different types of initial conditions 
(drivers) lead to rapid development of a turbulent state 
close to equipartition. During the following stage we 
model the free decay of these turbulence configurations. 
We study the development during this free decay stage 
and analyze characteristic parameters of the slowly vary- 
ing power spectrum of the turbulence. We assume that 
the properties of the turbulence configurations modeled 
in our numerical experiments are similar to those of cos- 
mological primordial turbulence. 

Our simulations allow us to estimate the spectral in- 
dices of the large-scale distribution of the kinetic and 
magnetic field energies in developed turbulence. It seems 



that the final configuration of primordial magnetic field 
depends on the initial conditions that must be deter- 
mined by the phase transition physics. On the other 
hand, the kinetic energy density of the developed state 
does not retain any information about the initial condi- 
tions and hence is not sensitive to the details of the phase 
transition. The spectral index of the kinetic energy at 
large scales can be well approximated as 2 (white noise 
spectrum), while the magnetic energy density spectral in- 
dex ranges between 2 and 4 depending on the initial con- 
ditions. Similar spectral indices for the large-scale mag- 
netic field are well established for laboratory turbulence, 
see Ref. |55|. Note that simulations in a finite periodic 



domain may suffer from the fact that originally discon- 
nected and causally independent regions come into causal 
contact within one box turnover time, which, based on 
the scale of the domain, is (vo^i) -1 - For the case shown 
in Fig. [U this time is about 15 ti, but the effective time 
can be even shorter owing to the effects of acoustic and 
Alfven waves. 

Our numerical results allow us to estimate typical time- 
scales for decaying free turbulence. It seems that cos- 
mological turbulence decays slightly slower then classi- 
cal grid turbulence in laboratory experiments. Although 
we took a rather small damping wavenumber, we can- 
not expect the establishment of a Kolmogorov-like spec- 
trum at small scales. We see a fast decorrelation of tur- 
bulence at small scales. We note that the phenomeno- 
logical approaches developed in Refs. 0, [28[ imply a 
fast Kraichnan-like decorrelation of turbulence. Conse- 
quently, only the large-scale or peak-scale magnetic field 
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results may have cosmological significance and contribute 
observable signatures, such as gravitational wave gener- 
ation and/or CMB anisotropy production. 

We have considered the case of non- helical turbulence. 
The presence of weak initial helicity can significantly 
change the development of turbulence. This is because 
magnetic helicity is a conserved quantity in the limit of 
large magnetic Reynolds numbers and can inverse cas- 
cade to larger scales [H, which could be cosmo- 
logically significant. This process has been confirmed 
through numerous simulations [25|, [26|, |3l|, E3] • 

Summarizing, we find that the generation of a mag- 
netic field at phase transition scales will lead to the de- 
velopment of turbulent motions, and, in the case of the 
electroweak phase transition, this turbulence has an ob- 
servable signature in the form of a gravitational wave 
signal. 
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